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ABSTRACT 

Future cosmic microwave background (CMB) polarisation experiments aim to measure an 
unprecedentedly small signal - the primordial gravity wave component of the polarisation 
field B-mode. To achieve this, they will analyse huge datasets, involving years worth of time- 
ordered data (TOD) from massively multi-detector focal planes. This creates the need for fast 
and precise methods to complement the M-L approach in analysis pipelines. In this paper, we 
investigate fast map-making methods as applied to long duration, massively multi-detector, 
ground-based experiments, in the context of the search for B-modes. We focus on two alter- 
native map-making approaches: destriping and TOD filtering, comparing their performance 
on simulated multi-detector polarisation data. We have written an optimised, parallel destrip- 
ing code, the DEStriping CARTographer (descart), that is generalised for massive focal 
planes, including the potential effect of cross-correlated TOD 1// noise. We also determine 
the scaling of computing time for destriping as applied to a simulated full-season data-set for 
a realistic experiment. We find that destriping can out-perform filtering in estimating both the 
large-scale E and B-mode angular power spectra. In particular, filtering can produce signifi- 
cant spurious B-mode power via EB mixing. Whilst this can be removed, it contributes to the 
variance of B-mode bandpower estimates at scales near the primordial B-mode peak. For the 
experimental configuration we simulate, this has an effect on the possible detection signifi- 
cance for primordial B-modes. Destriping is a viable alternative fast method to the full M-L 
approach that does not cause the problems associated with filtering, and is flexible enough to 
fit into both M-L and Monte-Carlo pseudo-C^ pipelines. 
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1 INTRODUCTION 

The temperature anisotropies of the cosmic microwave background 
radiation (CMB), and recently the polarisation anisotropies, have 
been used to constrain a set of cosmological parameter s to establish 
a standard model of the universe ISmoot et al. 1992.iHananv et all 



200d.lMelchiorri et al.l200(]llKovac e t al.'2002. Spergel et al.|2003l 
Hinshavy et al.ll2009llBrown et alj |2009b. Chia ng et alj|2009f) . The 



emphasis has now moved from revealing the contents and structure 
of the universe to explaining the unknown physics that generated 
them. Underpinning the standard model is the paradigm of infla- 
tion, which hypothesises a period of super-luminal expansion at 
very early times and identifies Gaussian, scale-free quantum fluc- 
tuations as the seeds of large-scale structure. Inflation models are 



predicted to uniquely generate tensor perturbations to the metric in 
the form of a stochastic gravitational-wave background that gener- 
ates a divergence-free "B-mode" polarisatio n pattern in the CMB 
that is not generated by scalar perturbations i Seliak & Zaldarria^ 



ll997l ; lKamionkowski. Kosowskv & Stebbinsll997h . The amplitude 
of the B-mode signal is inflation model dependent but is generally 
predicted to be orders of magnitude smaller than the temperature 
anisotropy. 

In order to detect such a subtle signal, we must build experi- 
ments of much greater sensitivity than has been possible so far. Due 
to the physical limits on the sensitivity of individual detectors, the 
next generation of CMB polarisation experiments aim to achieve 
this by using very large arrays of detectors. For example, the Q/U 
Imaging ExperimenT (QUIET) is operating in its first phase with 
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91 detector horns, each of which produces 4 independent measure- 
ments (2 of Q and 2 of U), and is planned to move to a 1500 horn ar- 
rangement in its second phase. Other examples include the balloon- 
born SPIDER, which will ope rate wit h 1024 detectors i n its largest 
band ( lMacTavi"sli et al. |2008|) . ebex jOxlev et al.l2005h which will 
use 14 40 detectors ove r three bands and the ground based POLAR- 
BEAR jLee et al.ll2008l) . which will use a total of 1274 detectors. 
The trend is the same for high-resolution temperature anisotropy 
experi ments, such as APEX- SZ, which is operating with 330 de- 
tectors jHalverson et alj2009l) . the Atacama Cosmology Telescope 
(ACT), which has begun operating w ith 1024 detectors per band 
( iKosowskvl l2003l : iFowler et alj l2007h , and the Sout h Pole Tele- 
scope, which has 465 in its largest band (Plagg e et al.ll2009h . 

To detect B-modes using the new generation of experiments, 
tight control of experimental systematics will be required. One of 
these systematics is correlated noise. Detectors display long term 
noise drifts that are characterised by a 1// power spectrum, re- 
sulting in correlated noise in the time-streams. In addition to this, 
ground-based experiments see an even larger 1// noise sourced 
from atmospheric fluctuations, wh ich may have a polarised ele- 
ment jHananv & Rosenkranzll2003l) . These noise sources can also 
be correlated between detector time-streams, for example through 
detector read-out electronics and due to the spatial correlation of 
atmospheric fluctuations projecting onto the focal plane. 

It is requisite of data analysis methods to characterise and re- 
duce the effects of noise and systematic errors. The pipeline in- 
volves a number of steps of radical data compression. The raw out- 
put from the telescope detectors is time-ordered data (TOD), a mas- 
sive dataset that contains signal varying with the telecope pointing 
information. TOD is compressed into a dataset of manageable size 
and in the natural form for a sky varying signal - a sky map. From 
the sky map, we estimate the set of angular power spectra, including 
the decomposition of the polarisation field into E and B-mode am- 
plitudes (accounting for spatial systematics such as foregrounds), 
from which we estimate cosmological parameters. The best place 
to remove correlated noise and time- variant systematic errors is in 
the map-making step. 

The two principle approaches to map-making for CMB ex- 
periments are the maximum-likelihood (M-L) meth ods (Tegmark' 
1997al . iDore et alj 1200 iL IStompor et all l200l Ide^G asperis et al, 
200a), which produce optimal maps slowly, and T OD filtering. 



which is used with Monte -Carlo pseudo-Cf methods l lHivon et al.l 
I2OO2I ISzaDudietal]|200lh and is very fast but necessarily sub- 
optimal. 

The M-L algorithms produce optimal maps in that they min- 
imise the noise whilst completely preserving the signal. The whole 
pipeline from time-ordered data to power-spectrum can be op- 
timal and maxim um-likelihood, including TOD noise e stimation 
(Ferreira & Jaffe 2000) and power-spectrum estimation l lTegmarlJ 
[l997h , ,Borrill,1999i) . the latter of which uses the analytic pixel- 
noise covariance matrix. This approach has been used very suc- 
cessfully and can be extended to include information on systemat- 
icsjStompor et al. 2002). Contemporary M-L codes, suc h as R OMA 
jde Gasperis et alj|2005h and MAPCUMBA jDore etal]|200lh . do 
not require explicit evaluation of the pixel-noise covariance matrix 
or its inverse (known as the "weight matrix"), and have been con- 
vincingly demonstrated as viable options for single-detector Planck 
data analysis including polarisation ( Ashdown et al. 2009, and ref- 
erences therein). A M-L code that includes detector correlations, 
called SANEPIC, has recently been developed and s uccessfully ap- 
plied to the BLAST dataset jPatanchon et al.ll2008h . a short dura- 
tion balloon experiment involving hundreds of detectors. 



As future data-sets become very large, the M-L algorithm be- 
comes increasingly difficult to implement. The method may be ten- 
able for some massively multi-detector experiments, but will re- 
quire massive and expensive computing platforms and very long 
computation times. The later stages of CMB data analysis re- 
quire information on noise power and correlation and the propa- 
gation of systematics in the maps that is effectively supplied by 
passing hundreds of Monte-Carlo simulations throu gh the map- 
making pipeline (eg: ^ Hivon et al.^ 2002^: MacTavish et al. 1 1200 8l : 
iBrown et alj|2009al : iTakahashi et aLlHoOT) . The M-L algorithm is 
not suited for such tasks, creating the need for fast methods that 
can be used on medium-sized computing platforms and used over 
and over again. Such methods could be used to complement the 
M-L algorithm in analysis pipelines, or even to replace it if their 
performances become close to optimal. 

The common alternative to M-L map-making is to apply a 
high-pass pre-whitening filter to the TOD. The noise part of fil- 
tered TOD is uncorrelated, so a map with uncorrelated noise can be 
returned by averaging the filtered TOD corresponding to each pixel 
(a naive map). This fast method is very well suited to Monte-Carlo 
simulations, but c ritically distorts th e signal part of the TOD. The 
MASTER method l lHivon et al]|2002l) mitigates this effect by evalu- 
ating a filter transfer function from Monte-Carlo signal-only reali- 
sations that can be deconvolved in multipole space. This approach 
is non-optimal and introduces extra variance into the power spec- 
trum estimates, especially at low-i" where the primordial B-mode 
signal is expected. 

The destriping me thod is being developed as a "third way" 
to CMB map-making teurigana et al.' 'l99T, 'Delabrouill^ 1 1998 , 



2005 



Mainoetal. 1999, Keihanen et al. 2004, Keihanen et al 
Kurki-Suonio et alj l2009l iKeihanen et alj 120091) . Destriping pre- 



whitens the TOD by approximating the low frequency noise part 
as a series of offset functions, which are then subtracted from the 
TOD. Long baseline noise drifts from 1// noise are subtracted 
without filtering the signal part of the TOD. The method is faster 
than M-L map-making and is tuneable in that the offset func- 
tion length can be varied, resulting in fast and dirty maps or slow 
and near optimal maps. Most destriping codes exploit the great- 
circle scanning strategy of satellite experiments to subtract long 
offset f unctions that are averages of the TOD in one re-pointing 
(see eg: [ Ashdown et al. 2007a for a discussion). The MADAM code 
jKeihanen et al. 2005.) . was the first to incorporate noise informa- 
tion, through a prior on the offset function amplitude distribution 
that permits the use of short offset amplitudes. A variant of destrip- 
ing that does not use noise information has recently been applied to 
total-intensity data from ACT to produce Sunyaev-Zeldovich (SZ) 
effect maps of galaxy clusters, using destriping baseli nes to model 
the in ter-detector common-mode atmospheric noise ( iHincks et alj 

I2OO9I) . 

In a previous paper jSutton et ai] I2OO9I) , we showed that 
destriping with short baselines using noise information is near- 
optimal when applied to non-circular scanning strategies, such as 
the cross-linked, azimuth only scans of ground-based experiments. 
In this paper, we extend the destriping algorithm to include in- 
formation on noise correlations between detectors and investigate 
the performance of our new code, the DEStriping CARTographer 
(DESCART), as applied to multi-detector TOD simulations. For 
these simulations, the Q and U signals are not modulated to the 
high-frequency white noise part of the spectrum, as is the case 
for half-wave plate modulation experiments, and are contaminated 
with 1// noise. We compare it to the MASTER method with fil- 
tering, the established fast map-making method, to determine the 
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strengths and limitations of the fihering approach for reahstic fu- 
ture data-sets. 

Given the importance of computation time for large data- 
sets, we also investigate the algorithmic scaling of DESCART as 
applied to a large full-season multidetector data-set. We concen- 
trate on polarisation measurements for B-mode experiments. The 
DESCART algorithm is also applicable to multi-detector SZ experi- 
ments, which we plan to investigate in a future paper. 

This paper is organised as follows: in S|2] we review map- 
making formalism of filtering and destriping in the context of mult- 
detector polarisation data; in !j3] we directly compare the filtering 
and destriping pipelines through simulations for a range of noise 
scenarios, comparing errors in timestreams, maps and E and B- 
mode angular power spectra; and in SQ] we apply DESCART to a 
much larger simulation of a massively multi-detector experiment, 
focussing on algorithmic scalings and the partitioning of the analy- 
sis. 



2 DESTRIPING AND FILTERING 
2.1 The map-making problem 

The map-making problem requires a solution for a sky map, Xp, 
given detector time-streams dt, pointing information and possibly 
noise information. The time-ordered data (TOD) is modelled as 



dt = AtpXp + fit 



(1) 



where Atp is a pointing matrix that contains all of the relevant 
pointing information, t indexes time and p indexes sky pixel. The 
noise time-stream nt is commonly assumed to be Gaussian and 
piece-wise stationary, satisfying 



(nt) = 
{ntrit') = Cn 



(2) 
(3) 



In general, nt is correlated noise, such that ^ is not diagonal. If 
the noise is stationary, its covariance ([3) is well described as sym- 
metric Toeplitz, or as a function of separation Cjvtf = f{t ~ t'). 
If the instrumental beams are symmetric, the pointing matrix for a 
temperature measurement Atp is composed of Is and Os and con- 
tains only one non-zero element per time row indicating which sky 
pixel the beam centre falls upon. 

This basic framework can be extended to multiple detectors 
and to measurements of polarisation. In the experimental configu- 
ration we simulate in this paper, each detector "pixel" makes direct 
measurements of Q and U Stokes parameters in its own frame of 
reference. The rotation from sky to detector frame of reference is 
subsumed into the pointing matrix 



RiOt) = 



Ftp — R{dt)Atp, 

cos 26t sin 29t 
— sin 26t cos 29t 



(4) 
(5) 



For simplicity in notation, we define the pointing matrix for 
detector-Q time-streams as Ptp = (cos 29 1, sin 26 1) Atp and for 
detector-U time-streams as Ptp = {— sin 29 1, cos 29 1) Atp. From 
this point we drop the distinction between them and consider P^p 
as the appropriate pointing matrix for the time-stream indexed by 
k. The map on which the polarisation pointing matrices operate 
contains maps of the Q and U Stokes parameters, concatenated end- 
to-end X = {x^ ,Xp). 

For multiple time-streams, we stack the polarisation pointing 



matrices, data and noise time-streams end-to-end to produce the 
multi-detector analogue to the data model, 
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(7) 



Given the data model, we can construct estimators for the sky 
map X. The M-L solution to the map-making problem is 



f = (P'^Cj 



'c-'^d 



(8) 



where P^ is the transpose of P and has the effect of summing 
TOD into a map, rather than projecting a map onto TOD. In the 
limit that we consider the TOD noise to be white (uncorrelated), 
the solution becomes very simple, corresponding to naive binning 
and averaging. 



(pTp-j-lpTJ- 



(9) 



2.2 The destriping solution 

The maximum-lik elihood method can be approximated by the de- 
striping algorithm jSutton et alj2009llAshdown et al.l2009l and ref- 
erences therein). The noise vector n is modelled as being composed 
of two components: uncorrelated (or white) noise; and a series of 
discrete offset functions that represent all of the correlated noise. 

With representing the random Gaussian realisation of 
white noise, the new data model for a single detector is 



dt=Pt 



tp^p 



-nY + Ft. 



(10) 



where the matrix Fta maps a set of basis functions, with amplitudes 
aa, onto the time-stream. 

For multi-detector systems, we generalise the basis func- 
tion amplitude vectors by concatenating them into a vector a = 
(a„ , , . . . , ) , and defining a block-diagonal multi-detector off- 
set pointing matrix 



/ F}^ 


V '■■ 





F?a 



The simplest choice of basis function is a constant offset 

Ftc = <! ^ * ^ 
° otherwise. 



(11) 



(12) 



where Aq is a chunk of the time-stream with length Ad. With this 
definition, Ftada becomes a vector of constant offsets, with am- 
plitudes a, approximating the correlated noise. More complicated 
basis functions c an be chosen, for example Fourier series or Legen- 
dre polynomials jKeihanen et al. ' 2004), though we focus on using 
short uniform baselines. 

The amplitudes a will be Gaussian random numbers satisfying 





Ca 



(13) 
(14) 



With a as a second set of parameters, the posterior for both the map 
and the offset amplitudes is 



P(x,d\d) oc P(d|f, d)P(x)P{d) 



(15) 
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where we have used the fact that the CMB map and the offset am- 
pMtudes are completely independent. We do not wish to assume any 
prior for the CMB mapQ, so we consider P{x) to be constant. 

The next step is to decide what prior information on the ampli- 
tudes to include through P{a). If we have an estimate of the power 
spectrum of the correlated noise, which the offset functions are de- 
signed to approximate, we can include prior information through 
the offset covariance matrix Ca. The probability distribution for 
the offsets is Gaussian 

P(a) = |Car'/%xp(^- ia^CaV^, (16) 

where Ca is the covariance of the offset amplitude estimates (ig- 
noring factors of 2n). If we have prior knowledge of Cn, we can 
build a prior Ca through 

Ca = (F'^F)~^F'^CnF(F'^F)"\ (17) 

The joint posterior is a product of two Gaussians: the likeli- 
hood of the data and the prior for the amplitude estimates. We can 
maximise the posterior by minimising the function /, where 

/ = -21n|P(f,a|d)| 

= (d-Fa-Pf)^Cw~^(d-Fa-Pf) 

+a^Ca^a. (18) 

Solving df /dx = 0, we obtain an expression for the map 

f = (P^Cw'P)~'P'^Cw'(d-Fa). (19) 

This result can be understood as follows: The correlated noise Fa is 
subtracted from time-stream d, leaving a time-stream composed of 
only signal and white noise, which is naively averaged into a map 
(the maximum-likelihood map for a time-stream with uncorrelated 
noise is the pixel average). 

To determine the offsets, we begin by substituting il9i into 
l |18t and gathering all the terms involving the pointing matrix P 
into a new operator, Z, 

/ = (d - Fa)^Z'^Cw Z(d - Fa) + a^Ca ^a, (20) 

where 

Z = I-P(P^Cw'P)"'P^Cw' (21) 

and I is the identity matrix. Z is a signal cleaning operator that 
estimates the naive map from the timestream (including both signal 
and any noise that is indistinguishable from signal) and removes it. 

Solving df /da — 0, we obtain an estimator for the offset 
amplitudes 

(F^Cw'ZF + Ca')3 = F'^Cw'Z(l (22) 

This is an inverse problem like that of the maximum likelihood 
algorithm. The system can be solved quickly, using the precondi- 
tioned conjugate gradients method, at a fraction of the processing 
time required to solve for the full maximum-likelihood map. An 
effective preconditioner K, can be constructed and applied easily 
from the offset prior covariance matrix 

K = (F'^CwF + Ca') (23) 

^ It is standard practice to assume priors on cosmological parameters, but 
the sky map is our first estimator of signal and should be free of map priors, 
which themselves have a very non-linear dependence on the cosmological 
parameter set. 



Destriping solves the map-making problem by making an ap- 
proximate model for the noise. In the limit that the offset func- 
tions perfectly model the correlated noise, the map solution is the 
maximum-likelihood map. 



2.2.1 Destriping correlated time-streams 

For multi-detector systems with noise uncorrelated between time- 
streams, the noise covariance matrix, Cn, will be non-zero only in 
diagonal detector-detector blocks (which themselves will be sym- 
metric Toeplitz). Cross-correlated noise will have non-zero off- 
diagonal detector-detector blocks in its Cn which will lead to non- 
zero diagonal block in Ca through equation l |17t . 

The inversion of the Ca is achieved by using some simplify- 
ing assumptions. To a good approximation, the individual detec- 
tor blocks can be considered to be circulant and can be inverted 
individually very easily. A circulant matrix becomes diagonal in 
Fourier space and each Fourier mode uj, which is independent, can 
be inverted as a scalar. In the presence of non-zero off-diagonal de- 
tector blocks, each mode becomes an independent N'^etector ma- 
trix, which we invert using Cholesky decomposition. This is the 
same approach a s is used in inverting th e multi-detector noise ma- 
trix in SANEPICi FPatanchon et alj ( l2008h . This entails the inversion 
of rioffscts matrices of dimension ridet each, once, before the itera- 
tions begin, an O (rioffscts n^^t) operation. 

The resultant matrix ((7^'° still has diagonal detector 

kk' blocks in Fourier space, so the matrix multiplication Ca "^a are 
conducted in Fourier space to minimise the operation count 

Wh{uj) = ^ Affcfc/(a;)?;fc/(a;). (24) 
fc' 

Here the vectors w and v are multi-detector vectors of offsets, and 
the matrix M is either Ca ""^ for the matrix application or K for the 
preconditioning step. 



2.3 The filtering approacii 

Filtering approaches to map-making apply a time-domain filter to 
TOD in order to suppress noise or systematics in the time-stream, 
before naively binning the now uncorrelated timestreams. For a fil- 
ter h(t — t'), the map estimate is 

f = (P'^P)"^P'^/i*d! (25) 

In this paper, we consider the effect of aggressive filters in- 
tended to suppress correlated noise rather than filters specific to 
experimental systematics. One choice of filter to suppress corre- 
lated noise is the prewhitening filter h = awC-^ , where aw 
is the white noise standard deviation. The covariance of the TOD 
filtered in this way is diagonal {a^^Sij). Another common choice 
is the overwhitening filter, h — a^C^ filter, which further sup- 
presses low frequency noise. 

The filter convolution is applied in the Fourier domain, where 
Cn is diagonal. This requires a Fourier transform pair to be ap- 
plied to the timestream, both of which have an operational scal- 
ing of 0{Nt log2 A'sogment ), where Nt is the number of TOD and 
A^'scgmont is the length of the segments over which the noise is sta- 
tionary and correlated (for example the length of a constant eleva- 
tion scan). 
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2.4 Covariance of the estimated maps 

Both of the pipelines used in this paper construct the final map us- 
ing a naive binning step. The noise covariance of maps constructed 
using naive binning is 

C = (P'^P)~^P'^Cn°°P(P'^P)~\ (26) 

where Cj^^ is the covariance of the TOD. For the filtering 
pipeline, C^f^ = J e'^'^ hCNh^ duj, which for the prewhitening 
filter is diagonal and produces a diagonal pixel covariance. 

In the destriping pipeline, two sets of parameters are solved 
for; the map and the offset amplitudes. The Fisher matrix for the 
full parameter set is 



d^f/dx^ d^f/dxda 

d'-^f/dadx d^f/da^ 

P^C^^P F^C^^P 

p'^c; 



(27) 



where / is the posterior in l |18t . Inverting this matrix by partition 
gives the noise covariance of the destriped maps 



We minimize the memory usage of the code with a juggling of 
the data. When all the timestreams that might be correlated (those 
within a single scan) are loaded they are immediately used in two 
ways: added to an accumulated naive map, and reduced into base- 
line offsets. They are then discarded. The application of the prior 
and likelihood matrix operations during the PCG requires only a 
single allocated one-detector timestream at a time (per process). 
Finally, when the best-fit offsets have been computed, they are con- 
verted directly into a map (which represents the correlated noise in 
the original naive map). Deducting this from the first map yields 
the output map. At no point do we allocate full timestreams. With 
this modification DESCART can process 900 one-hour timestreams 
per 2-GB RAM process for typical baseline lengths. This increases 
markedly if multiple time-streams are associated with the same 
horn. 

As currently constituted the code is not much slower than gen- 
erating naive maps for larger runs, since simply loading the data 
takes a significant proportion of the run time. 



C = (P'^Cw'P- (28) 

(P^Cw'F)(F^Cw'F + Ca')-^(F^Cw'P))"\ 

which can be solved fo r low-resolution maps of large angular scales 
jKeskitalo et alj2009l) . 

This solution assumes that the destriping data model is correct 
- that the correlated noise is perfectly described by the offset func- 
tions. Whilst this is an acceptable approximation for short offset 
function lengths, it breaks down for long baseline functions. In this 
regime, unmodelled correlated noise contributes significantly to the 
map covariance. 

The quantity that destriping tries to measure is the average of 
the noise in a segment of the TOD, often referred to as the "refer- 
ence baseline", defined by r = {F'^ F)^^ F'^ dn where d„ is noise- 
only TOD. The covariance described by equation ([28} is from the 
error between the baseline amplitudes estimated by l l22t and these 
reference baselines. If this error is uncorrelated with the unmod- 
elled correlated noise within a TOD segment, as it is for our sim- 
ulations (see ^3.2\ , we can account for the unmodelled noise by 
adding a correction term C*^""""" to the covariance matrix l |28l l, 

(-<corr ^ (pTp)-lpTcdesp^pTp^-l ^29) 

Cn " = Cn + FCaF'^ 

-CnF(F'^F)"^F'^ + F(F'^F)"^F'^Cn, (30) 

where Ca is the offset prior and Cn is the TOD noise covari- 
ance. For the symmetric Toeplitz Cn and Ca we consider, Cf^^ 
is also symmetric Toeplitz. The implementation of i29\ requires an 
0{NtNcorr) summation, where Neon- is the correlation length of 
the reference-offset subtracted timestream. Any correlation in this 
timestream is due to unmodelled 1// and has a short correlation 
length. 

2.5 The code 

The DESCART code in which we implemented the above algorithm 
is written in Fortran and uses MPI for communication; this ensures 
scalability to the largest available systems. The communication op- 
erations occur when generating and distributing the collective naive 
map for the current iteration, and when correlating horns in a scan. 



3 COMPARISON OF THE PIPELINES 

In this section, we analyse the comparative performance of the 
filtering and destriping pipelines as applied to simulations of a 
constant elevation scanning multi-detector telescope. We begin by 
analysing time-stream and map domain errors, the latter for which 
we use the root-mean-square (rms) residual as the comparative 
statistic. We proceed to analysing spherical harmonic domain sys- 
tematics and errors, for which we need to estimate the E and B 
mode angul ar power spectra. F or this task we use the Monte-Carlo 
approach of'Hivon et al.' ('2002)), generalized for polarization power 
spectra as in Smith ( 2006) such that estimators do not mix E and B 
modes. E to B mixing is mitigated by applying an apodizing win- 
dow function to the estimated maps. The window function we use 
is a symmetric circular cosine apodisation around the map cen- 
tre. The apodization is the same as that used in our previous pa- 
per (Sutton et al. 2009), to which we refer the reader for details. 
The method requires estimation of the power spectral noise bias 
(iV/^^) through Monte-Carlo (MC) noise only simulations and es- 
timation of the filter transfer function Fi through signal only sim- 
ulations. Finally, we analyse the statistical significance of the B- 
mode detection provided by the two pipelines. 

Figure [T] shows the range of filters investigated in the filtering 

— 1/2 

pipeline. These comprise prewhitening (Cp^ ) and overwhiten- 
ing (Cj,^""") filters, for low frequency noise suppression, and a scan 
frequency filter, used for the suppression of scan-synchronous sys- 
tematics. 

The destriping results in this section make use of inter- 
timestream noise correlation information. The benefit of includ- 
ing these correlations is dependent on the experimental set-up and 
for the polarised dataset we simulated, the noise reduction when 
including the correlation was small. For other systems, including 
this information may be important. When applying the algorithm 
to simulated data in ^ correlation information was ignored as the 
experiment has negligible inter-detector correlated 1//. The choice 
of destriping length is informed by the noise fk'- as we are in- 
corporating prior noise information, a choice of destriping length 
Ad < 0.1/ fk is used to to approximate the optimal M-L map as 
closely as possible. 
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Figure 1. Filters used in the filtering pipeline. The soUd curves are the 
— 1/2 

prewhitening filters Cj^ and the dashed curves are overwhitening fil- 
ters Cj^"*". Colour denotes the noise case, with black denoting fi^ = 5-mHz 
a = 2, blue f)^ = 50-mHz o = 1, green /j, = 200-mHz a = 1, and 
red fk = 200-mHz a = 2. The dotted curve is the scan frequency filter, 
defined for the scanning frequency and its first harmonic with a bandwidth 
of 0.02-Hz. 
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Figure 2. Plan of the simulated focal plane arrangement. Each horn pro- 
duces a Q and a U time-stream in the detector's frame of reference. 



3.1 Simulated Data-sets 

Time-ordered data was simulated for a variety of noise cases, 
shown in Table [T] using a one-day scanning strategy comprising 
four constant elevation scan sets. This scan is a typical observa- 
tion of a single non-galactic field from an observatory in the Chan- 
jnantor Scientific Reserve in Chile. At a sampling frequency of 
fa = 100-Hz, the total integration time for the scans were 6 hours 
and 38 minutes. We simulated an instrument with a Gaussian sym- 
metric beam with FWHM= 12' and a focal plane consisting of 19 
horns arranged in a hex pattern, each producing 2 time-streams. We 
show this arrangement in Figure[2l 

For each noise case, we generated 100 signal-plus-noise and 
100 independent noise-only realisations. Each realisation consisted 



Table 1. Noise parameters used in the simulations. The spectral index of 
unity describes the correlated noise from detectors, whilst the spectral index 
of two desciibes atmospheric noise. The ranges of knee frequency chosen 
are representative of frequencies experienced by low frequency (/ < 100- 
GHz) experiments. The and a of noise case 6 are representative of 
QUIET Q-band data. The values of aw here correspond to one year's of 
observing time with 1000, 91 and 19 horns, where cr^r has been scaled so 
as to produce the white noise amplitude in the final map expected for these 
dataset sizes. For all cases, we simulate weakly correlated inter-detector 
f~°' noise, with correlation coefficient p = 0.2. 



Parameter 
Symbol 


White noise RMS 
crw 


Knee frequency 
fk 


Spectral index 

a 


Case 1 


130-^1 A' 


200-mHz 


1 


Case 2 


59.3-p,K 


200-mHz 


1 


Case 3 


\1.9-p,K 


200-mHz 


1 


Case 4 


\1.9-p.K 


200-mHz 


2 


Case 5 


\1.9-p.K 


50-mHz 


1 


Case 6 


\1.9-pK 


5-mHz 


2 



of 38 time-streams, with Gaussian cross-correlated 1/f noise, 
Gaussian white noise, and an independent CM B signal realisation 
produced using SYNFAST ( Gorski et ai]|200f ) from the input Ces 
used in our previous paper jSutton et al.l2009l . for which the tensor 
to scalar ratio r = 0.1). 

Cross-correlated noise was generated by statistically corre- 
lating the 1/f streams. For each time-stream, a Fourier-space 
Gaussian random number vector 4>d{f) was generated, where 
d indexes detectors. We built a correlation matrix A^^i — 



Pdd'<^d<Td'\/ fk.dfk,d'' where ad and fk,d are the noise standard 
deviation and knee frequency for stream d respectively and pdd' 
is the correlation coefficient for streams d and d' . Correlated 1/f 
noise streams were produced by applying the Cholesky decompo- 
sition of this matrix to the set of ipd mode-by-mode 



nd{f) = P^'\f)Y.Ldd'<t^d'{f), 



(31) 



where P{f) = (1//)" and A = LL^. 

To get time-domain noise stream realisations with the desired 
auto and cross power spectra, each nd{f) stream was Fourier trans- 
formed and added to an independent Gaussian white noise realisa- 
tion with standard deviation ad- The white noise level itself was 
determined assuming a fiducial noise-equivalent Q (NEQ) of 248- 

3.2 Time-streams and maps 

Both pipelines approach map-making by suppressing noise in the 
time domain before naive binning. In the upper panel of Figure [3] 
we show the results of using the methods on mean time-domain 
noise power The time-stream power spectrum is composed of sig- 
nal (red curve) and noise (black curve). Filtering suppresses both of 
these components, enforcing the flatness of the output TOD power 
(green curve). Destriping only affects the noise component of the 
TOD, suppressing its power at frequencies smaller than the destrip- 
ing offset function length (blue curve), in this case at 1-Hz. 

The noise in the destriped TOD comes from two components: 
unmodelled noise and destriping error. Unmodelled noise includes 
all noise sources at frequencies higher than the destriping length 
and some power beneath it, which decreases exponentially at fre- 
quencies beneath the destriping length, behaviour well described 
by the application of a transfer function to the noise spectrum, as 
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Figure 3. Upper panel: Power spectra of TOD for 1// noise (black curve), 
filtered noise (green curve) and destriped noise (blue curve - the destriping 
length corresponds to / = 1-Hz). Also plotted is an arbitrarily normalised 
power spectrum of signal only TOD (red curve). Lower panel: Power spec- 
trum of unmodelled noise (black curve), power spectrum of destriping error 
F(!*(red curve) and the sum of these spectra. 



described bv lKurki-Suonio et alj200^ . Defining the reference base- 
lines as r = (F'''F)~^r^d„ (ie: the average of the noise within 
a baseline) the unmodelled noise is d„ — Fr. The destriping error, 
defined by F(a — r), describes the error in estimation of the base- 
line amplitude. This error is a complicated quantity, determined by 
the pointing of the telescope. 

The lower panel of Figure[3] shows the power spectrum of the 
unmodelled noise (black curve), the power spectrum of the destrip- 
ing error (red curve) and the sum of these power spectra (green 
curve). The destriping error is the dominant component of the noise 
in the destriped TOD at low frequencies. At high frequencies, the 
unmodelled noise dominates, whilst the destriping error shows a 
series of troughs at multiples of the destriping length caused by the 
window function of the offset functions. 

The sum of the power spectra is virtually identical to the power 
spectrum of the destriped noise TOD. For our simulations, this val- 
idates the assumption, made in i]2.4| that the destriping error is un- 
correlated with the unmodelled noise. Putting the unmodelled noise 
du ~ dn — Ff, and the destriping error e = a — r, the power spec- 



Table 2. RMS map residuals (in /iK). See Table[T]for a description of the 
noise cases. The filter considered is the prewhitening filter. The destriping 
length was = 1-s for all except the ff^ = 5-mHz simulations (noise case 
6), where it was Ad = 20-s. 



noise case 


descart Q 


descart U 


filtered 


filtered U 


1 


7.46 


7.41 


6.91 


6.86 


2 


3.40 


3.38 


3.21 


3.20 


3 


1.03 


1.02 


1.17 


1.19 


4 


1.01 


1.00 


1.18 


1.20 


5 


0.98 


0.97 


1.04 


1.04 


6 


0.94 


0.93 


0.99 


0.99 



trum of the destriped TOD is 



/ 



dte 



-ift 



{{du-Fef)= I dte"^'{{dl) + {{Fef) 

-{due^F'^)-{Fedl)), (32) 



which requires the cross-correlations 



'F^) - (Fedl) 



0. 



This facilitates the use of l |30| l to correct the destriped map covari- 
ance matrix as calculated by J29t . 

Figure |4] shows example Q and U maps from DESCART and 
their residual maps (defined as estimated map minus the map used 
to simulate the data). The maps from the filtering pipeline are 
not shown, as they appear identical by eye. (For well cross-linked 
strategies, 1// is a subtle systematic that appears most strongly in 
power spectra). The comparative power in the filtered residual maps 
is dependent on the signal to noise ratio. Filtering forcibly produces 
white noise in the output maps, optimally suppressing correlated 
noise. In noise dominated maps, this makes filtering more effective 
in minimising map residuals. 

We compare the residual maps, which give the best measure 
of errors on the sky, by comparing the residual root mean square 
(rms) over the realisation ensemble. Table [2] shows the rms resid- 
ual for both pipelines for all the parameters simulated. When the 
signal to noise ratio is small (19 and 91 horns), the residuals from 
filtering are smaller than those for destriping. For the higher sig- 
nal to noise simulations (1000 horns), destriping produces smaller 
residuals than filtering. By a quirk of the parameterization, 1//^ 
noise produces less noise at frequencies higher than fk than does 
1// noise, leading to larger residuals in noise case 3 compared to 
noise case 4. 

Filtering suppresses signal as well as noise and in high sig- 
nal to noise regimes, the signal distortion contributes appreciably 
to the residuals. Signal distortion contributes more to the residuals 
when the required filtering is more aggressive, when either the knee 
frequency or spectral index of the correlated noise is increased. Ex- 
ample maps of signal distortion are shown in Figure|5] The magni- 
tude of the signal distortion is not dependent on the integration time 
over a pixel, but rather on the level of cross-linking and polarisation 
angle dispersion in the pixel. 



3.3 B-mode systematics from filtering 

The systematic error effects of T OD filtering are wel l understood in 
the case of temperature data (eg: iHivon et"al]|2002h . the dominant 
effect being the convolution of the estimated power spectrum with 
a filter transfer function, Ff . This can be determined from MC 
signal-only realisations and can be deconvolved from the estimated 
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Figure 4. Example maps of Q and U Stokes pai'ameters in /iK from DESCARTfrom simulations with white noise for 91 horns (/fe 
The top row is the map estimate and the bottom row is the residual map. A 10° graticule is overlaid. 



200-mHzanda = 1.0). 



power spectrum. In the case of polarisation, E and B-mode filter 
transfer functions can also be measured through simulations. 

Our simulations show that B-modes estimated using the filter- 
ing algorithm are biased by a spurious B-mode component caused 
by E— ^B leakage. This leakage is absent when no filtering is ap- 
plied to the TOD, as is the case with the destriping algorithm. It is 
also absent when the input maps are contrived to have no E-modes. 

Figure [6] shows the spurious B-mode from signal-only TOD 
simulations with zero input B-mode. The left plot shows contrast- 
ing effects of TOD filtering: suppression of the observable B-mode 
power (by convolution with the filter transfer function Ff); and 



spurious B-mode power caused by distortions of the Q and U maps 
(such as those seen in Figure[5](. Both the B-mode suppression and 
the spurious B-modes are greater for more aggressive filters with 
higher knee frequencies and spectral indices. The spurious B-mode 
is non-negligible for all of the noise cases considered at the band- 
powers of interest for measuring the primordial B-mode (which 
peaks at f ~ 100), with an amplitude of the order of the observable 
signal for r = 0.1. 

The right hand plot of Figure |6] shows the spurious B-mode 
and signal suppression from overwhitening and scan-frequency fil- 
ters in comparison to the prewhitening filter. The overwhitening 
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Figure 5. Signal distortion in /iK from (prewliitening) filtering in signal-only maps. The left panel is Q and right panel is U. 




Figure 6. E— J-B leakage systematic induced by filtering. The left panel shows spurious B-mode power (error-bars) compared to the model B-mode power 
spectra (solid curves) suppressed by the prewhitening filter transfer functions. The least aggressive filter (red curves) suppresses the input B-mode power 
spectrum the least and also distorts the Q and U maps the least, resulting in less E— >B leakage and less spurious B-mode power The red plots are for filtering 
with /(j = 50-mHz and spectral index « = 1, the black plots are for = 200-mHz and a = 1, the blue plots are for fi^ = 5-mHz and a = 2, and the green 
plots are for = 200-mHz and a = 2. The right panel shows these effects for the overwhitening and scan frequency filters. Colour denotes the noise case 
(blue is case 4 and red is case 6), with the solid curves and plain error bars showing the prewhitening filter, dashed curve and diamond error bars showing the 
overwhitening filter, and the dotted black curve and star error bars showing the scan frequency filter. 



filters produce more spurious signal (and to a lesser degree, more 
signal suppression) than the equivalent prewhitening filter for the 
same noise case. The scan frequency filter produces significant spu- 
rious signal, particularly at large angular scales, where it is larger 
than the observable B-mode signal. 



The right hand plot of Figure |7] shows the biasing effect of 
the spurious B-modes on Ff as measured from MC signal-only 
simulations. Ff is the ratio of the Ces of the filtered signal only 



map (including signal bias) and a binned map of signal only TOD 



Ff 



rpBonly , 7-1 + 



(33) 



where Sf is the spurious B-mode signal that leads to an additive 
filter bias, F^, of the filter function. The true B-mode filter trans- 
fer function is shown in the right panel of Figure |7] (dotted curve, 
estimated from B-mode only signal simulations), compared to the 
additive filter bias (dashed curve, estimated from filtered Simula- 
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tions with zero input B-mode) and the resulting biased function 
(sohd curve). 

The power spectra can be debiased by extending the MASTER 
approach to include a spurious signal bias in the same way as it 
includes a noise bias. If Ki yy is the band-power cou pling kernel as 
defined for polarisation in [Hansen & Gorskil j2003l) . the unbiased 
power spectrum Ci can be returned by 

(d) = K~,]Py,{Ct - (iVf ^) - {5f ^)) (34) 

where Pu is the binning operator defined in lHivon et al. I( l2002h and 
A'^ and Si are noise and spurious signal biases calculated over sets 
of MC simulations. 

The E B leakage contributes more than just power . The 
variance of the leaked E-modes contributes to the variance of the 
debiased B-mode estimates (see ij3.4| >. 

The results presented in this section are limited to the single 
experimental set-up described in i ]3.1| We have found that the bias 
and its variance can vary with the details of the scan-strategy and 
focal plane arrangement, and consideration should be given to min- 
imising this effect when implementing filtering methods. 

3.4 Polarisation power spectra errors 

Using both pipelines, we have estimated unbiased mean E and B- 
mode signals for each of the noise cases in Table [T] In the filtering 
pipeline, the modified estimator l l34t was used to remove spurious 
power from E ^ B leakage . Figure [8] shows an example of the 
mean E and B-mode estimates over the ensemble of realisations for 
one of the parameter sets together with bandpower variances calcu- 
lated over the realisation ensemble. In this section, we evaluate the 
pipelines based on the bandpower variances of the debiased E and 
B-modes. 

TOD filtering decreases the sensitivity of the power spectra 
to signal, particularly at large angular scales, as described by the 
filter transfer functions, Fi, shown in Figure [7] for the different 
noise cases and filters. These must be deconvolved from the power 
spectra estimates to produce an unbiased estimator of the CMB. 
Part of the gain of filtering noise correlations before mapping is 
counteracted by this loss of sensitivity to signal, as the signal-to- 
noise in heavily filtered modes is larger than can be achieved by 
optimal methods. In our simulations, destriping consistently pro- 
duces smaller variance in the power spectra than filtering. Destriped 
power spectra do not require an Fi deconvolution, as the estimator 
is by construction unbiased with respect to the sky signal and pro- 
duces flat filter transfer fuctions. 

Figure |9] shows the ratio of the E-mode bandpower standard 
deviation from filtering to that from from destriping, presented as a 
percentage error increase, and its variation with the noise correla- 
tion parameters, and a, for both prewhitening and overwhitening 
filters. The error increase is larger for the more aggressive filters 
with higher knee frequency (or alternatively, filtering for simula- 
tions with more correlated noise). Likewise, the more aggressive 
overwhitening filter results in larger bandpower variances at large 
angular scales than the corresponding prewhitening filter for the 
same noise case. 

The error increase is present in all signal to noise regimes sim- 
ulated and shows only weak dependence on the signal to noise ratio 
in the maps. Figure [TO] shows the effect of varying signal to noise 
on the E-mode error increase (for noise cases 1-3, with fk = 200- 
mHz and a — 1). This contrasts with the map pixel variance, where 
filtering and destriping outperform one another in noise and signal 
dominated regimes respectively (see Table[2]l. As the filter transfer 
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Figure 9. E-mode error increase when using filtering instead of destrip- 
ing, for simulations with varying /j. and a. The solid curves ai'e for the 
prewhitening (Cj^'*'''^) filter and the dashed curves ai'e for the filter 

function is scalar, it is not dependent on the magnitude of the white 
noise, but on the noise correlation statistics. It does not depend on 
the ratio of the signal degradation in the map to the noise in the 
map and so the effect of white noise variation on the power spectra 
should be small. 

B-mode bandpower errors are affected by both the loss of sen- 
sitivity due to the TOD filter and the extra variance from E B 
leakage described in Section [331 Figure [TT] shows the B-mode er- 
ror increase for the different noise cases and filters. The increase 
is considerably higher than for the E-mode, amounting to > 50% 
error increase at large angular scales for all of the noise cases sim- 
ulated. Low frequency noise filtering has a more pronounced effect 
on the error increase for B-modes than for E-modes, resulting in 
larger error increases for higher knee frequencies, higher spectral 
indices (in contrast to the E-mode, as discussed above) and for the 
overwhitening filter. The disparity between this and the E-mode can 
be explained by the greater influence of the E ^ B leakage, for 
which variations in spectral index have a stronger effect. As noted 
in the previous section, the size of leakage, and thus the error in- 
crease reported here, can vary depending on the experimental setup. 

3.5 B-mode detection significance 

The primary goals of near-future CMB polarisation experiments 
are the detection and characterisation of the lensing component of 
the B-mode power spectrum, and the potential detection of the pri- 
mordial gravity wave component. To achieve this we will require 
high precision from our analysis methods and any extra error they 
contribute to the angular power spectrum will have an impact on 
these science goals. 

To constrain the effect of bandpower variance differences on 
B-mode detection, we investigate the B-mode detection signifi- 
cance achieved by the pipelines. Figure [12] shows the total B-mode 
detection significance per bandpower they return. 

The hump like shape of these curves is due to increasing band- 
power variance due to Gaussian white noise at small angular scales 
and due to decreasing signal power at large angular scales. Destrip- 
ing produces higher significance bandpower detections, with the 
improvement increasing when the correlation of the noise (through 
fk and q) is higher, or if the filter is more aggressive. 
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Figure 7. Effects of noise correlation parameters on E and B-mode filter transfer function. The left plot shows the E-mode function for the filters considered. 
The solid curves denote prewhitening filters and the dashed curves denote overwhitening filters, with the colour of both denoting the noise case. We also show 
the filter transfer function for the scan-frequency filter The right plot shows the bias effect of the spurous B-modes on the filter transfer function F^-^ for the 
case with /j. = 200mHz and a = 1. The dotted curve is the unbiased filter transfer function p^°"^y from B-mode only simulations, the solid curve is the 
biased B-mode transfer function and the dashed curve is the additive filter bias F^'^ . 
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Figure 8. Mean E-mode (left) and B-mode (right) estimates for simulations with /j. = 200mHz and a = 1.0, with white noise for 1000 horns. The black 
crosses are destriped estimates and the red crosses are filtered estimates. The solid lines are the input fiducial power spectra and the dashed fine in the B-mode 
plot is the lensing B-mode component of the input spectrum. The large-scale B-mode bandpower variances are considerably enhanced for filtering compared 
to destriping, in this case preventing detection of the primordial B-mode peak. 



A directly comparable statistic can be calculated by combin- 
ing the bandpowers into a single data point, measuring the total 
significance of the detection. For this we use the total signifi cance 
estimator defined in our previous paper l iSutton et al.l 120 09*), that 
produces the single data point by binning weighted to a fiducial 
model (7/"* - in this case the known input B-mode power spectrum 
- and weighted inversely by the the bandpower variance erf 



C = 



(35) 



Table [5] shows the B-mode detection significance, defined as 



[c / ^ {[c - {c)Y)]. 

We also show the significance of the detection of the primor- 
dial B-mode peak in the bandpower centered &t I = 100. Whilst 
our chosen input noise level has not supported a 2a detection for 



Table 3. Total significance of primordial plus lensing B-mode detection. 



fk 



a descart C, 



-1/2, 



5-mHz 


1.0 


13.44 


12.49 


12.12 


50-mHz 


1.0 


13.05 


11.56 




200-mHz 


1.0 


12.64 


10.69 




200-mHz 


2.0 


11.77 


9.98 


9.44 



the peak, it shows filtering can shave almost Icr off the weak detec- 
tion from destriping. 
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Figure 12. B-mode bandpowers detection significance (in multiples of the bandpower error a) for different noise correlation parameters, /(; and a, and signal 
to noise. Top left: fh = 50-mHz, a = 1.0 Top right: = 200-mHz, a = 1.0. Bottom left: fk = 200-mHz, a = 2.0. Bottom right: = 5-mHz, a = 2.0. 
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Figure 10. Effect of signal to noise on E-mode error increase. For simu- 
lations with white noise level for 19, 91 and 1000 detector horns, using 
fk = 200-mHz and a = 1. Only N"-'^/'^ filtering is considered. 



Table 4. Significance of B-mode inflationary peak detection. 



/j; a descart C'jv^^^ filter C^^^ filter 



5-mHz 


1.0 


1.87 


1.28 


1.18 


50-mHz 


1.0 


1.83 


1.20 




200-mHz 


1.0 


1.72 


1.05 




200-mHz 


2.0 


1.78 


0.87 


0.70 



4 APPLICATION TO A MASSIVELY MULTI-DETECTOR 
EXPERIMENT 

The speed of map-making methods will become vital in future 
CMB experiments, as data-set sizes increase. It is important there- 
fore to understand the algorithmic scalings of the methods available 
and whether destriping is competitive as a fast map-maker. Destrip- 
ing is linear with the number of pixels in the map, making it useful 
as a method for high resolution analysis, but its scaling with data- 
set size is the more relevant question for massively-multidetector 
small field experiments. 

In this section, we address this question by applying the 
DESCART code to simulated full-season data from a massively 
multi-detector experiment to determine the scalability of the 
method to much larger datasets than those in the previous section. 
In ij4.1l we determine the scaling of computing time and iteration 
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Figure 11. B-mode error increase when using filtering instead of destrip- 
ing, for simulations with varying and a. The sohd curves are for the 
prewliitening (Cj^^''^) filter and the dashed curves are for the Cf^^ fil- 
ter. The error increase is larger than for the E-mode due to the additional 
variance from E— >B mixing by TOD filtering. 



number with the dataset size. In i]4.2| we investigate the effect on 
the map pixel variance of partitioning the analysis by scan, a pro- 
cedure that may be necessary for very large datasets. 

The simulations were produced using the pointing information 
from the first full-season of QUIET Q-band operation. The sim- 
ulated dataset comprises of 202 constant-elevation scans (CESs) 
of a CMB sky patch with a total wall-clock observation time of 
704.6 hours at a sampling frequency of 50-Hz. Each CES includes 
timestreams from 68 detetors - produced by 17 hexagonally ar- 
ranged horns with four polarisation sensitive diodes per horn. This 
produced 4.7 * 10* hours of integration time, summed over all the 
detectors. The total memory requirement for the whole dataset was 
67-Gb, rising to 101-Gb including pointing information. 

1// noise was simulated to be uncorrelated between detectors 
and between CESs. We choose a noise of 200-mHz and an offset 
length of Ad = 2-s. These simulations were produced to determine 
the scaling of software and use typical, fiducial noise characteris- 
tics. As such, neither the noise characteristics nor the simulations 
they produced are similar to real QUIET data. 

4.1 Algorithmic scalings 

We have investigated the scaling of computing time with data-set 
size by applying DESCART to an increasingly larger set of CESs. 
In this scheme, all of the included data is used to estimate and sub- 
tract the naive map of the offset solutions - the operation ZFa from 
equation J22l l. Convergence of the PCG iterator is achieved when 
the offsets from all of the included time-streams satisfy the conver- 
gence critereon - that the error vector ||e|| = \\Ax — fe|| < t||6||, 
where || || denotes the 2-norm and the arbitrary stop tolerance r = 
10~®. Using all of the available data to estimate the offset vector a 
means that the size of the system to be solved becomes very large 
when realistic datasets are involved. 

However, the computation time of the algorithm scales lin- 
early with the quantity of input data. Figure [T3] shows the scaling 
of computing time with the number of CESs included in the analy- 
sis. The computing time does not include the initial data-load oper- 
ations, which will be system architecture dependent, and the time 



plotted is an average: the total computing time required to analyse 
all 202 CESs for the upper panel, and 40 CESs for the lower panel, 
divided by the number of partitions of that size required to analyse 
the whole data-set. In our analyses of the maps, the majority of the 
total computation time was taken by the data load. 

The scan sets shown in both panels of Figure [T3] were pro- 
cessed with DESCART on the 610 node cluster Titan, at the Univer- 
sity of Oslo0. The larger scan sets in the upper panel were pro- 
cessed using 64 cores on 8 nodes, whilst the smaller sets of the 
lower panel were processed using 16 cores on 2 nodes, so the pan- 
els are not directly comparable. 

The operation count per iteration is dominated by the signal 
subtraction operation l l21l l that includes two large linear 0{Nt) 
processes - map binning and map projection. These dominate the 
smaller 0{Na log A'^q) EFT operations associated with the prior 
Ca, as Na is two orders of magnitude smaller than Nt- The total 
computation time depends on the number of iterations taken by the 
PCG algorithm to return a sufficiently precise estimate of the off- 
set amplitudes. For the scanning strategies we consider, the number 
of iterations does not increase as more data is added: rather it de- 
creases. This behaviour is shown in Table[5] which shows the mean 
iteration number required per CES set. This can be understood by 
noting that increasing the number of CESs included in the analysis 
provides a better estimated naive map in the Z operator of | |2U . 
so the operation ZFa becomes closer to Fa (an effect enhanced 
by the lack of correlation in noise between CESs). With less noise 
wrongly identified as signal, the condition number of the destriping 
matrix in J22b decreases and so the number of iterations required to 
solve the simpler set of equations reduces. 

The reduction in the iteration number is most marked when 
CESs are considered in pairs instead of independently. This is due 
to the increased level of cross-lirtking in the combined scanning 
strategy of the two CESs, both in revisitations of the map pixels and 
in variation of the polarisation-sensitivity angles for the pixel. The 
vast majority of the reduction in iteration number in our simulations 
is achievable by including groups of 4 CESs. 

The above arguments can be applied to increases in the num- 
ber of horns. In the absence of inter-hom correlated noise, further 
horns are computationally identical to further CESs (other than in 
terms of cross-linking). We therefore argue that the scaling of the 
algorithm with horn numbers will also be linear. 

No inter-timestream cross-correlated noise was included in 
the simulations for this section, so the effect of including a cross- 
correlated prior as in i]2.2.1l was ignored. Including this effect 
results in an additional operation per iteration, in which the in- 
verse of the Fourier space inter-detector prior matrix is applied to 
the Fourier space data (equation I24t. which does not change the 
number of FFTs applied each iteration. The operation scales as 
0{ncorrNjNa), where Na, Nd and the numbers of de- 

striping offsets, detectors and correlated pairs for each CES. In the 
case that the number of detectors and correlated detector pairs are 
high, this operation could come to dominate the 0{Nt) operations 
discussed above. 

4.2 Analysis Partitioning 

The memory requirement for analysing complete data-sets be- 
comes untenable for massively-multidetector ground-based exper- 
iments, so partitioning the data for the destriping step is an attrac- 

^ http://hpc.uio.no/index.php/Main_Page 
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Figure 13. Linear scaling of computation time with number of CESs in- 
cluded in destriping analysis, with a dotted linear curve shown for refer- 
ence. Upper panel: Mean computing time per CES set for whole data-set 
(202 CESs) using 64 processors. Lower panel: Mean computing time per 
CES set for smaller dataset (40 CESs) using 16 processors. 



Table 5. Variation of mean iteration number with the number of CESs in- 
cluded in the analysis. For CES numbers between 40 and 202, the number 
of iterations remains constant. The mean iteration number is defined as the 
total number of PCG iterations required to map the whole dataset divided 
by the number of CES sets required to analyse it. 



#CESs 


# Iterations 


1 


27.25 


2 


17 


4 


16.3 


8 


15.4 


40 


14 



such that Zn = n, then the offsets become the reference offsets 
described in ^'i.2\ We can expect that the pixel noise in the destriped 
maps will reduce as the offset amplitude estimates better model the 
correlated noise in the timestreams. 

Destriped maps were produced for different destriping parti- 
tions, then combined to produce a final map that uses all of the 
available data (in this case, the 40 CES data-set). For a set of maps 
Xi, with Q/U white noise covariances Ci, the combined map is de- 
fined as f EE (E, Cj-i)-^ Ci-^xl. 

This ensures that the noise amplitude in the final maps is con- 
stant and any differences between them are solely due to changes 
in the error of the offset amplitude estimates. Figure [14] shows the 
variation of the RMS residual of Q and U in the final combined map 
with destriping partition size. Only pixels within the central region 
of the scan (the science field itself) were used to calculate the RMS 
residual. 

The effect of partition size on RMS residual amounts to a 
0.3% reduction between destriping with 40 CESs and destriping 
with single CESs only. Of this change, the majority of the improve- 
ment is achieved by considering pairs of CESs - as was the case 
with iteration numbers and for the same reasons. The residual RMS 
remains constant for CES numbers up to the full 202 CES dataset. 

A number of effects that naturally produce cross-linking 
within a single CES contribute to make this change small. In a 
separate analysis, we considered data from the central horn only, 
to remove the cross-linking effects of the 17 horn focal plane. The 
variation in RMS with CES number is qualitatively similar to that 
shown in Figure [74l except that the RMS reduction is much larger 
- amounting to a 3% reduction in residuals when all 40 CESs are 
used to estimates offset amplitudes rather than considering each 
CES separately. 

The effect of destriping analysis partition ing has been inves- 
tigate d recently for the Planck experiment by ( iKurki-Suonio et aU 
I2OO9I) , who found a significant variation in destriping pixel residual 
RMS with the length of scan considered. However, we note that the 
scanning strategies of the Planck and QUIET experiments have lit- 
tle in common - the Planck scanning strategy used revisited the vast 
majority of pixels only after a 6 month spin period, with the destrip- 
ing relying on the cross-linking in a small number of polar pixels to 
determine amplitudes for offsets the length of a single great circle 
scan. The scanning strategy considered here is heavily cross-linked 
for most pixels on small time-scales. 

Our simulations suggest that destriping partitioning is an ac- 
ceptable analysis technique when the partition consists of a few 
CESs. A caveat for this result is that it is based solely on the 
effects of 1// noise. Instrumental systematics, the simulation of 
which is beyond the scope of this paper, can be better constrained 
with larger data partitions. An example of this is scan-synchronous 
noise, which can be modelled by the offsets when large numbers of 
CESs with uncorrelated scan-synchonous noise are included, and 
which significantly increases the number of iterations required to 
solve for the offset amplitudes. 



five analysis option. Destriped maps can be combined with a noise- 
weighted average - an operation that is mathematically identical, 
other than through the accuracy of the offset estimates, to making 
a map from the entire destriped dataset. 

However, the accuracy of the offset amplitude estimation can 
be increased by including more data in the solution of equation [22] 
As argued in i]4.1| the signal removal operator Z, which removes 
both signal and any noise that looks like signal, removes less noise 
if more data is considered. If it were to remove no noise at all. 



4.3 Resources Summary 

A year of data consisting of 600 hours of data with 70 detectors 
sampled at 50Hz can be analyzed with this code on 64 typical mod- 
em processor cores in approximately 30 minutes. This will scale 
approximately linearly with increasing time and number of detec- 
tors. 

This is fast enough that destriping can be included as part of a 
Monte- Carlo pipeline, replacing map-making with heavy filtering 
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Figure 14. Variation in RMS residual for the Q and U Stokes pai'ameter 
maps with destriping partition size. The maps were combined such that the 
same quantity of data was used for each map. The variation seen is therefore 
entirely due to offset estimation error, as opposed to white noise level in the 
maps. 



with only a small increase in computing time. It is also suitable as 
an alternative to maximum-likelihood mapmaking for performing 
repeated auxiliary tasks such as null tests in cases where running 
the full pipeline repeatedly is unfeasible. 

The code and method can therefore be applied to the largest 
presently existing CMB polarization data sets to make high resolu- 
tion maps. 



5 DISCUSSION 

The maximum-likelihood (M-L) map-making method has long 
been used in CMB data analysis pipelines to reduce sizeable time- 
ordered data-sets (TOD) into optimal maps, which are optimal in 
the sense that the noise in the map is minimised without the loss 
of information - the sky signal in the maps is not distorted. How- 
ever, applying the method is moving toward being untenable for fu- 
ture long duration, massively multi-detector CMB experiments, es- 
tablishing the need for faster, approximate, map-making methods. 
Such methods will be crucial to simulating and removing experi- 
mental systematics with enough precision to search for the small 
primordial gravity wave signa l in the CMB polarisa tion field. 

Significant study (see eg: lAshdown et alj2009l and their refer- 
ences) has been conducted into data-analysis for space based exper- 
iments, such as the Planck experiment. We have built upon this by 
evaluating fast map-making methods for massively multi-detector 
ground base d experiments, of w hich a large number are in devel- 
opment (see lBrown et alj|2009al for a review). We have developed 
DESCART, an optimised, parallel destriping code, and applied it 
to simulations of time-ordered data (TOD) from such an experi- 
ment. We destripe using short baselines with a noise prior - a mode 
of operation that h as been shown to produce near optimal maps 
dSutton et alj|2009l) . 

For large future datasets, the fastest map-making method un- 
der consideration is TOD filtering. Our comparison of the filtering 
and destriping approaches shows that, for the highly cross-linked 
scanning strategies we simulate, TOD filtering underperforms in 
power spectrum errors when compared to destriping. This result 
is due to two effects: the suppression of signal sensitivity by the 



noise filter, which decreases the signal to noise ratio in the lower 
I bins, despite removing the effects of correlated noise; and the 
introduction of E ^ B mixing by the TOD filter, which can be 
characterised and removed by TOD simulations, but which signif- 
icantly contributes to the variance of the large angular scales. Of 
these effects, the latter (only present in the B-mode spectrum) is 
much more dominant, typically doubling the bandpower variance 
at the expected inflationary B-mode peak at ^ ~ 100, although the 
variance increase is dependent on how aggressive the TOD filter is. 

In our simulations, we find that bias from filtering can strongly 
affect the detection of the primordial B-mode peak (see Table |4]l. 
This is true for the one realisation of possible scan strategies, patch 
shapes and receiver array arrangements presented here; the full 
variation of the effect with the parameters of the scan and exper- 
iment is beyond the scope of this paper, and we expect that it can 
be ameliorated in many situations. It is also unclear how the leak- 
age effect will scale for longer scanning-strategies and larger focal 
planes. 

DESCART's computing time depends on the length of the 
dataset and the condition number of the destriping matrix. Each it- 
eration of DESCART is dominated by operations that scale linearly 
with the size of the dataset - the P and F operators and their trans- 
poses. The number of iterations is sensitive to the condition of the 
matrix, which is generally improved by the addition of more heav- 
ily cross-linked data, such that the number of iterations required to 
solve for offset amplitudes reduces slightly. For the simulations we 
considered, the computing time was dominated by the initial data- 
loading operation. This suggests that destriping can be potentially 
as fast a method as filtering. 

The effect of an extended focal plane is to add sufficient cross- 
linking to a single scan that destriping offsets are well constrained 
by data from that scan alone. The gain from including multiple 
scans to estimate offset amplitudes is a reduction of order 0.3% in 
map residuals, a gain that is almost entirely returned by combining 
a small number of scans. 

The search for B-modes will likely be dependent on how well 
experimental and physical systematics can be constrained and miti- 
gated. Many of these systematics require aggressive filtering of the 
time-stream data prior to map-making. However, destriping has the 
potential to model some of these systematic effects in addition to 
correlated time-stream noise. For ground-based, constant-elevation 
scans, the most important systematic is scan-synchronous signal, 
caused for example by ground pick-up in the side-lobes of the ex- 
perimental beam and time-varying atmospheric noise. The pres- 
ence of scan-synchronous signal tends to add considerably to the 
number of PCG iterations required to solve for offsets, as offsets 
have a limited capacity for modelling it. An alternative to filtering 
these modes is to solve for additional destriping offsets mapped 
by azimuth rather than by time. Such an approach is equivalent 
to the estimation and removal of ground s ignal recently applied 
to the QUaD data by iBrown et alj ( l2009bh . For total power ex- 
periments, the offsets can be mapped onto all the detectors within 
the horn, much as atmospheric common-mode offsets are mapped 
onto timestreams f rom multiple adjacent horns in the method of 
iHincks et ai]( l2009h . 

As part of the generalization of destriping to multiple detec- 
tors, we have included the possible effect of inter-detector corre- 
lated 1 // noise in the offset prior. The use of this information in 
the prior results in a marginal improvement in destriping perfor- 
mance, but for the experimental set-up we simulate, the effect is 
small. We note, however, that in other experimental situations, the 
inclusion of the correlation information could be important. 
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We have investigated the use of correlation information to 
constrain cross-correlated 1//, but have ignored the possibility of 
cross-correlated white noise. Correlated white noise tends to af- 
fect detectors within the same horn, in which case it is optimally 
accounted for by including the correlations in the white noise co- 
variance matrix, which we have considered to be uncorrelated here. 
The detectors in the same horn observe the same sky pixels, so this 
extra information does not lead to replacing the naive map evalu- 
ated in each destriping iteration with a more difficult solution, as the 
naive mapping operation is still diagonal in pixel-space. We plan to 
address experimental systematics, including correlated white noise 
and scan-synchronous signal, in a future paper. 
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